function [MSE,LS,CRPS] = crossval(vadata,vamdl)

numTCH = numel(vadata); 
VAdist = vamdl.VAdistEstimates;
numCOMP = vamdl.NumMixtureComponents;

ERR = cell([numTCH 1]);
LS = cell([numTCH 1]);
CRPS = cell([numTCH 1]);

parfor j = 1:numTCH

    n = numel(vadata(j).Y);
    
    ERR_j = zeros([n 1]);
    LS_j = zeros([n 1]);
    CRPS_j = zeros([n 1]);
    
    for i = 1:n

        vad = vadata(j);
        vad.Y(i) = [];
        vad.Z(i,:) = [];
        vad.ZF(i,:) = [];
        vad.sid(i) = [];
        [vahat,vapost] = vaposterior(VAdist,vad);

        Y = vadata(j).Y(i);
        Z = vadata(j).Z(i,:);

        ERR_j(i) = Y - Z*vahat;

        Y_mean = Z*vapost.mean;
        Y_var = arrayfun(@(c) Z*vapost.cov(:,:,c)*Z' + VAdist.residvar,1:numCOMP,'unif',1);
        Y_sd = sqrt(Y_var);

        LS_j(i) = log(normpdf(Y,Y_mean,Y_sd))*vapost.pr';
        
        z = (Y - Y_mean)./Y_sd;
        crps = Y_sd.*((1/sqrt(pi)) - 2*normpdf(z) - z.*(2*normcdf(z) - 1));
        CRPS_j(i) = crps*vapost.pr';
        
    end

    ERR{j} = ERR_j;
    LS{j} = LS_j;
    CRPS{j} = CRPS_j;

end

MSE = mean(cat(1,ERR{:}).^2);
LS = mean(cat(1,LS{:}));
CRPS = mean(cat(1,CRPS{:}));

   